An energy landscape model for glass-forming liquids in three dimensions 
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We present a three-dimensional lattice-gas model with trivial thermodynamics, but nontrivial 
dynamics. The model is characterized by each particle having its own random energy landscape. The 
equilibrium dynamics of the model were investigated by continuous time Monte Carlo simulations 
at two different densities at several temperatures. At high densities and low temperatures the model 
captures the important characteristics of viscous liquid dynamics. We thus observe non-exponential 
relaxation in the self part of the density auto-correlation function, and fragility plots of the self- 
diffusion constant and relaxation times show non-Arrhenius behavior. 

PACS numbers: 64.70.Pf 
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I. INTRODUCTION 

A mechanical system of N spherically symmetric par- 
ticles is completely characterized by its so-called energy 
landscape, the graph of the potential energy function 
t/(ri, ttv) in 3N + 1 dimensions. As suggested by 
Goldstein in his pioneering 1969 paper U, the energy 
landscape is particularly useful for elucidating the dy- 
namics of highly viscous liquids. This is because viscous 
liquid dynamics are dominated by jumps over barriers 
much larger than fc^T; most time is spent on vibrations 
around local energy minima of the landscape. However, 
it was only after the work of Stillinger and Weber in the 
1980's and the enormous growth in use of computer 
simulations in the 1990's that the energy landscape be- 
came a dominant paradigm in the study of viscous liq- 
uids 13, H, |E 0, IE IE EE El ■ For recent reviews see, e.g., 

Elil. 

It is difficult to imagine a complex high-dimensional 
landscape, but an obvious idea is to assume that there 
is an element of randomness in the landscape. In this 
philosophy one follows Wolynes, who argued that some 
phenomena occurring in a specific complex system are 
typical of those that occur in most systems chosen ran- 
domly out of an ensemble of possible systems . 

A possible disordered landscape consists of a high- 
dimensional lattice with random, uncorrelated energies 
chosen, e.g., according to a Gaussian, with nearest- 
neighbor Metropolis dynamics. This model, which has 
trivial thermodynamics, has been shown to reproduce a 
number of observed properties of viscous liquids, and the 
low-temperature dynamics of the model are understood 
to be dominated by site percolation However, 
the model does not have a meaningful thermodynamic 
limit; if the distribution of energies is chosen such that 
the mean energy is extensive, the relaxation times are not 
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intensive. The problem, of course, is that a single nearest 
neighbor jump on the lattice changes the energy by an ex- 
tensive amount, effectively corresponding to a complete 
rearrangement of all molecules. Another problem is how 
dimensionality is reflected in the energy landscape. Many 
condensed matter systems behave differently in two and 
three dimensions. If this applies also for viscous liquids, 
it must somehow be reflected in the landscape. 

The question we consider here is: Is it possible to con- 
struct a sensible 'generic' random landscape model? Such 
a model should obey the following requirements: 

1. It should have a well-defined thermodynamic limit, 
i.e., extensive average energy and intensive relax- 
ation times. 

2. It should reflect the dimensionality of space. 

3. All sites should be statistically equivalent, thereby 
ensuring translational invariance on the average. 



II. THE MODEL 

The energy landscape is attractive because it abstracts 
from three dimensions. Nevertheless, we would like to 
suggest that the simplest way to have a model obeying 
the requirements listed above is to return to three dimen- 
sions: 

Consider a lattice gas in three dimensions. If ran- 
dom energies are assigned to the lattice sites, the sys- 
tem is described by Fermi statistics. This corresponds 
to particles in an external random potential, thus with 
no translational invariance and only the trivial self exclu- 
sion particle-particle interaction. A simple modification 
turns this model into a highly nontrivial model, namely 
to assume that each particle has its own energy land- 
scape, i.e., the energy of the system is given by (where 
Ti is the lattice position of particle i, and S is the Dirac 
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(5- function): 



(1) 



The first term is tlie energetic interaction; for each 
lattice site, r, and each particle, the energy, ei(r), is 
chosen randomly from a probability distribution, p(e). 
ei(r) and are generally different random numbers 

if i 7^ J, i.e., particles have different energy landscapes. 
The second term in Eq. Q is the self exclusion particle- 
particle interaction; no more than one particle is allowed 
at each lattice site. 

Now all lattice sites are statistically equivalent. The 
model allows calculation of pressure and chemical po- 
tential, and has extensive thermodynamics and intensive 
relaxation times. In particular we expect that at high 
densities, p = N/V, there will be a jamming effect slow- 
ing down the dynamics considerably. 

When comparing the simulation results of this model 
to results from molecular dynamics (MD) simulations, 
one should keep in mind that the model does not include 
the high frequency vibrations associated with "cage- 
rattling" . The dynamics that are modeled here are the so 
called "inherent dynamics" [^, i.e., the result of map- 
ping the true dynamics onto a series of inherent struc- 
tures (local minima in the 3N — 1 dimensional energy 
landscape) . 

We make the simplest possible choice for the probabil- 
ity distribution p(e): the Box distribution \p[t) = 1,0 < 
e < 1] . In this case the mean system energy per particle 
is easily found to be (/3 = l/fc^T): 



CB) ^ 1 
N (3 
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(2) 



exp(/3) - 1 

At low temperatures, we thus get (E) = NksT. 

III. SIMULATION DETAILS 

The model was simulated on a three dimensional 
L X L X L cubic lattice using the N-fold way kinetic 
Monte Carlo Method 0, with continuous time. We 
use Metropolis transition rates with local Monte Carlo 
moves; if a particle jumping to a nearest-neighbor site 
brings the system from state i to state j, the associated 
transition rate is given by: 



T{i ^ j) = min[ro. To exp{-(3{Ej - E,))] 



(3) 



Our length and time units are defined by setting the lat- 
tice unit a = 1, and the fastest transition rate Fq = 1. 

Since each particle has its own energy landscape, the 
number of different site energies are given by x = 
p X L^. Storing these numbers in memory would put a 
severe constraint on how large systems we could simu- 
late. Instead we utilize the 'ran4' random number gen- 
erator 13 in the following way: each particle is assigned 
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ph = l- p 


/3- values 


L 


0.992 


8 X 10"^ 


0,2,4,..., 16 


10 


0.999 


1 X 10"^ 


0,1, 2,. ..,13 


20 



TABLE I: Parameters used in simulations, p = N/L'^ 
the density of holes (unoccupied sites). /3 = l/fcsT 



ph IS 



a 'particle-seed', and when needed this is used together 
with the appropriate site-index as input to 'ran4', which 
performs a series of bit operations to produce a uniform 
random deviate in the range 0.0 to 1.0. 

The simulations are carried out with two sets of pe- 
riodic boundary conditions, one for each term in eq. Q 
We denote by L the lattice length associated with the 
particle-particle interactions. The usage of 'ran4' makes 
it possible to use much larger energy landscapes: the lat- 
tice site index used as input to 'ran4' is a 32-bit integer, 
and for the energy landscapes we can therefore use a side- 
length of approximately 1000 (we use L times an integer), 
i.e. for all practical purposes (at least in this paper) each 
particle has an energy landscape that is infinite. 

As mentioned above, we impose locality on the Monte 
Carlo moves; particles can only jump to (vacant) nearest- 
neighbor sites. Relaxing this requirement (letting a 
Monte Carlo move consist of a random particle inter- 
change with a random particle/hole) gives us an efficient 
way to equilibrate the model: For all the state points in- 
vestigated here the characteristic time for equilibration 
was found to be less than 5 time units. Equilibration 
runs were done for a time period of 1000 time units. 



IV. SIMULATION RESULTS 

Two densities were simulated, each at a range of (3- 
values, see table ^ Our simulations agree with the ana- 
lytical expression for the mean energy (equation |2J|. Re- 
ported results are averages over 8 independent simula- 
tions (8 different energy landscapes), and error-bars are 
estimated from fluctuations between these 8 simulations. 



A. Mean-square displacement 

Fig. n shows the mean-square displacement, (r^(i)), 
in a log- log plot. These results look similar to what is 
found in MD simulations of viscous liquids (see eg. [20j): 
At long times the dynamics is diffusive {{r'^{t)) oc t), and 
this diffusive regime is preceded by a plateau that devel- 
ops as the system is cooled. In MD simulations this devel- 
oping plateau is attributed to " cage rattling" : particles 
vibrating in a cage consisting of the nearest neighbors. 
In this regime particles move considerably less than the 
inter-particle distance. This "MD scenario" is obviously 
not what happens in this model; as discussed above the 
vibrations on length scales shorter than the inter-particle 
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FIG. 1: Mean-square displacement, (r-^(t)), for p;, = 1 x 10 ^ 
(upper panel) and ph ~ 8 x 10~^ (lower panel). See table Q 
for further details. Data points are connected by straight 
lines. Error-bars indicate 95% confidence interval in {r^{t)) 
estimated from fluctuations between 8 independent samples 
(uncertainties only discernible at short times). 



distance are not included in the model. Here the develop- 
ing of a plateau means that after a particle has jumped to 
a nearest-neighbor site, the probability for jumping back 
to where it came from is (on average) larger than the 
probability for jumping to a new lattice site. This leads 
to a slowing down of the dynamics compared to diffu- 
sion dynamics. Only when particles have jumped several 
times is this correlation between jumps lost whereupon 
the dynamics become diffusive. At short times {t < 1) 
a regime with (r^(t)^ oc t is seen - this is simply a con- 
sequence of the time scale being so short that particles 
never jump more than once. 

In Fig. we report the diffusion coefficients extracted 
from Fig. ^ For reference we show here also results for 
the p — limit, i.e., simulations with non-interacting 
particles (in this limit the model is obviously not a good 
model of a liquid). In the p = limit we find Arrhenius 
behavior [D = Dq exp{~(3AE)], as expected from perco- 
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FIG. 2: Diffusion coefficients extracted from the mean square 
displacements in Fig. Q For reference the diffusion coeffi- 
cient in the p = limit (non-interacting particles) is included. 
Lower panel: Same data as upper panel. Y-axis scaled by 
Do = D(J3 — 0). X-axis scaled by /3g to make the data col- 
lapse at (3//3g = 1. At the density p^ = 1 x 10"'^ we define 
f3g — 13. For ph ~ 8 X 10~^ and p^ = 1 empirical scaling was 
used to find jSg — 13.95 and Pg = 30.9 respectively. For p = 
a straight line was fitted to the data. For the high densities 
data points are connected with straight lines. 



lation arguments [21|. In contrast, the higher densities 
show distinctive non-Arrhenius behavior; the model ex- 
hibits "fragile" behavior. To facilitate comparison with 
Angell's fragility plot [22|, we show in Fig. the dif- 
fusion coefficients scaled in the following way: the y-axis 
is scaled with the diffusion coefficient at infinite temper- 
ature, Dq = D{f3 — 0) (which scales with ph H^), and 
the X-axis is scaled with an "inverse glass temperature", 
Pg, which is here defined by the scaled diffusion coef- 
ficients being identical at f3/Pg = 1 (and (3g = 13 for 
Ph = 1 X 10""^). The degree of fragility is observed to 



FIG. 3: Apparent activation energies calculated from eq. |3 
Squares: ph = 1 x 10~^. Triangles: ph = 8 x 10~^. Circles: 
Ph 1 (i-e., p = limit, non-interacting particles). Data 
points are connected with straight lines. 



increase with increasing density (decreasing ph). 

Fig. 13 shows the apparent activation energy, obtained 
by regarding AE in the Arrhenius expression as being 
temperature dependent: 



AE = 



(4) 



In the p = limit AE as expected approaches the the- 
oretical value, Ec ~ 0.31 24J. At high densities, AE 
keeps increasing above this value, reflecting the non- 
Arrhenius behavior. There is an indication (particularly 
for pfi — 8 X 10~^) that AE might be levehng of to a 
constant, indicating that there might be a crossover from 
non- Arrhenius to Arrhenius behavior, as seen e.g. in sim- 
ulations of viscous silica ■ The observed indication of 
a crossover to Arrhenius behavior might be related to 
"hitting the bottom" of the energy landscape, but this 
point deserves further investigations. 



B. Density-density correlation (self part) 




In Fig. 21 we show the self part of the density auto- 
correlation, Gs{0,t) [2^, i.e., the probability that a par- 
ticle at time t is at the same site as it was at time t — 0. 
Dashed lines are fits to stretched exponentials: 

fit) = cxpi-it/rr) (5) 

The fits are not perfect, but they capture the main char- 
acteristics of the data. The fitting parameters are shown 
in Fig. El As for the diffusion constant (Fig. [JJ, the 
relaxation time r exhibits non-Arrhenius behavior with 
an indication of a crossover to Arrhenius behavior at 



FIG. 4: Upper and middle panel: Self part of the density- 
density correlation, Gs{0,t) for ph = 1 x 10~^ and ph = 8 x 
10^^ respectively. Full lines are straight lines connecting data 
points. Error-bars indicate 95% confidence interval in Gs(0, t) 
estimated from fluctuations between 8 independent samples. 
Dashed lines are fits to stretched exponentials, exp(— (t/r)^). 
Fits were done for Gs(0, t) < 0.8. Lower panel: 1 — Gs(0, t) 
for ph^8x 10"^. 
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FIG. 5: Fitting parameters from fitting stretched exponen- 
tials to the return probability, Ga{0,t) (fig. 0J. Upper panel: 
relaxation time r. Lower panel: stretching exponent 7. Data 
points are connected by straight lines. 



the lowest temperatures. The stretching exponent 7 de- 
creases with /3, indicating an increasing degree of non- 
exponential relaxation. Except for the lowest tempera- 
ture at each density, the stretching exponent 7 seems to 
level off to a constant close to 0.5. A constant stretch- 
ing exponent indicates time-temperature superposition 
(TTS), i.e., that the shape of the relaxation function is 
independent of temperature. We note that this behav- 
ior is consistent with experiments indicating that TTS 
is correlated to 7 = 0.5 Here we find at the very 

lowest temperatures an indication that the stretching ex- 
ponent starts to increase again, which might be related to 
the apparent cross-over from Arrhenius to non-Arrhenius 
behavior discussed earlier. Simulations at lower temper- 
atures are needed to investigate this question further. 

Comparing Fig. ^to Fig. ^ and Fig. one might 
ask: "why does a plateau develop in (r^(t)) and not in 
Gs{0,t)7" . The answer is, that there is indeed a plateau 
developing in Gs{0, t) - this can be seen in Fig. where 
we have plotted 1 — Gs{0, t) (i.e. the fraction of particles 
that are contributing to (r^(t))) in a log-log plot. In fact. 




FIG. 6: Displacement vectors, at a time where (r^) « 2, for 
/3 = and 14 respectively (p — 0.992, see table Hll for details). 
The initial position of particles that moved to a new lattice 
site during the time interval is indicated by a filled circle, and 
the displacement vectors are shown as straight lines. 



at short times where no particle jumps more than once 
(r^(t)) = 1 — Gs(0, i). At the lowest temperature in Fig. 
g|; this relation holds to within 10% for t < 100. 

We note that in the p = limit at low temperatures 
Gs{0,t) looks quantitatively different from what is seen 
in Fig. ^ There is a strong initial relaxation (related to 
jumps with AE < 0), a stronger stretching (7 w 0.3), 
and a final pronounced power law regime that starts at 
Gs(0,t) « 0.1. This limit (which we again stress is not 
a good model for a liquid) is investigated in a separate 
paper p3 |. 

C. Dynamical heterogeneity 

It is well established that viscous liquids contain dy- 
namical heterogeneities, i.e., if subsets of particles are 
defined by their dynamical properties, these tend to be 
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/3 




t 


Gs{0,t) 





2.02 


6.2 X 10^ 


18.6% 


14 


2.06 


2.0 X lO'' 


54.1% 



TABLE II: Parameters describing the two different sets of 
displacement vectors in Fig. |S| Note: Averages are here only 
over particles, not ensemble/time averages. 



correlated in time and/or space [23, y^- Figure |B| indi- 
cates in a qualitative way that the model exhibits dy- 
namical heterogeneity to an increasing degree as temper- 
ature is lowered. In Fig. |^ we show the displacement 
of particles at a time where <^r^(t)) w 2 for = and 
/3 = 14 {pii = 8 X 10"'^). It is evident from the figure 
that the fraction of particles contributing to the mean 
square displacement (i.e. 1 — Gs{0,t)) is smaller at the 
low temperature (see also table Hljl , and that the positions 
of contributing particles are correlated in space. 



diffusion coefficients has been proposed. The first results 
from simulations of the model have been presented. At 
high densities the model exhibits the two non's character- 
izing viscous liquids, non-exponential relaxation and non- 
Arrhenius temperature dependence of relaxation times 
and diffusion coefficients. The fragility increases with 
density. Furthermore, indications of a number of inter- 
esting features was found; i) Non-Arrhenius to Arrhenius 
transition, ii) Time-temperature superposition, iii) Dy- 
namical heterogeneities. 
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